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ABSTRACT 

We present results from self-consistent numerical simulations of cosmic structure for- 
mation with a multi-frequency radiative transfer scheme and non-equilibrium molec- 
ular chemistry of 13 primordial species (e~, H, H+, H~, He, He+, He"*"+, H2, D, 
D+, HD, HeH+), performed by using the simulation code GADGET. We describe our 
implementation and show tests for ionized sphere expansion in a static and dynamic 
density field around a central radiative source, and for cosmological abundance evo- 
lution coupled with the cosmic microwave background radiation. As a demonstrative 
application of radiative feedback on molecular gas, we run also cosmological simula- 
tions of early structure formation in a ~ 1 Mpc size box. Our tests agree well with 
analytical and numerical expectations. Consistently with other works, we find that 
ionization fronts from central sources can boost H2 fractions in shock-compressed gas. 
The tight dependence on H2 lead to a corresponding boost of HD fractions, as well. 
We see a strong lowering of the the typical molecular abundances up to several orders 
of magnitudes which partially hinders further gas collapse of pristine neutral gas, and 
clearly suggests the need of re-ionized gas or metal cooling for the formation of the 
following generation of structures. 

Key words: cosmology: theory - structure formation 



■ 1 INTRODUCTION 



The preset understanding of cosmic structure formation re- 
lies on the observations of a Universe expanding at a rate 
of Ho — 70km/s/Mpc and whose energy budget is largely 
dominated by a form of unknown 'dark' energy or cosmolog- 
ical constant. A, that contributes ~ 70% to the total cosmic 
energy content. The residual matter contribution is roughly 
~ 30%, but only a very small fraction of ~ 4% consists of or- 
dinary baryonic matter, while the rest is unknown cold (i.e. 
non-relativistic) 'dark matter' (DM). More precisely, recent 
determinations suggest no,m ~ 0.272, fio.A = 0.728, and 
f^o.b = 0.044 (Komatsu et al. 2011). In this framework, also 
called ACDM model, cosmological structures can grow from 
gravitational instability (Jeans 1902) of primordial matter 
fluctuations, probably originated during the primordial in- 
flationary epoch. These early perturbations represent the 
seeds which would develop into present-day galaxies and 
stars (Schwarzschild & Spitzer 1953) by gas cooling and con- 
densation (Spitzer 1962), and they will affect the surround- 
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ing environment through a number of mechanical, chemical 
and radiative processes commonly known as feedback effects 
(on this topic, see the extensive review by Ciardi & Ferrara 
2005). 

Quantitatively speaking, linear perturbation analyses are 
usually performed to study the initial phases of gravita- 
tional collapse, where a Gaussian density distribution for 
the primordial matter ffuctuations is assumed. The linear ex- 
pansion of the continuity, Euler, and energy equations can 
also be extended with higher-order corrections (e.g. Tseli- 
akhovich & Hirata 2010; Maio et al. 2011; Stacy et al. 2011; 
Greif et al. 2011) or non-Gaussian deviations (e.g. Grinstein 
& Wise 1986; Koyama et al. 1999; Komatsu et al. 2002; 
Grossi et al. 2007; Desjacques et al. 2009; Maio & lannuzzi 
2011; Maio 2011; Maio & Khochfar 2012), but to study non- 
linear regimes and feedback mechanisms it is essential to 
perform numerical integrations and use N-body/hydro sim- 
ulations. Indeed, to capture early gas collapse it is needed 
not only to follow gravity and hydrodynamics, but also its 
full chemistry evolution and molecule formation. Since in the 
cosmic medium hydrogen (H) is the most abundant species 
with a cosmological mass fraction of Xh — 0.76 (correspond- 
ing to 0.93 in number fraction), its contribution in gas 
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cooling is likely to play a very relevant role, together with 
helium (He). However, H and He coUisional processes are 
able to cool the medium to ~ 10* K via resonant line transi- 
tions, but they are not capable to bring the gas temperature 
further down. At such low temperatures, thermal collisions 
are not able to excite the electrons to higher levels due to the 
large energy gaps (from a few to some tens of eV) of H and 
He atomic configurations. Saslaw & Zipoy (1967) proposed 
that gas cooling and fragmentation could be continued below 
~ 10" K by H2 cooUng, down to 10^-10^ K. Later, Lepp & 
ShuU (1984) suggested that the existence of primordial deu- 
terium (D) could determine HD formation and consequent 
cooling even below ~ 10^ K, down to several ~ 10 K. During 
the last decades, these problems have been tackled by col- 
lecting full reaction networks (Shapiro & Kang 1987a; Puy 
et al. 1993; Galli & Palla 1998; Hui & Gnedin 1997a; Abel 
et al. 1997; Uehara & Inutsuka 2000; Nakamura & Umemura 
2002; Omukai et al. 2005; Glover & Abel 2008; Omukai et al. 
2010; Glover et al. 2010) and by running high-resolution 
chemistry cosmological simulations, both in the standard 
ACDM model (e.g. Abel et al. 2002; Bromm et al. 2002; 
Yoshida et al. 2003, 2007) and in dark-energy cosmological 
models (Maio et al. 2006). Further on, the effects of metal 
cooling have been investigated in numerical simulations by 
joining molecular chemistry evolution with metal pollution 
and low-temperature fine structure transitions from, e.g., 
C, O, Si, Fe (e.g. Maio et al. 2007). These studies clearly 
show the strong implications of metals on the cosmological 
chemical evolution (chemical feedback) and how the rapidity 
of early enrichment from first star formation episodes (sec 
also Tornatore et al. 2007) overcomes gas molecular cool- 
ing (Maio et al. 2010, 2011), marking the transition from 
the primordial, pristine star formation regime - i.e. the so- 
called population HI (popIII) regime - to the more standard 
population H-I (popH-I) regime (Maio et al. 2010). This 
transition is often parameterized in terms of a minimum 
critical metallicity Zcrit ~ 1O~'*Z0, but given our ignorance 
about early dust production and detailed atomic and molec- 
ular data, Zcrit has large uncertainties: its expected value 
varies between ~ 10"® Z© (Schneider et al. 2003) - 10"^ Z© 
(Omukai et al. 2010) and ~ IQ-^Z© (Bromm & Loeb 2003). 
There are also recent arguments on the possible existence of 
a critical dust-to-gas ratio (Schneider et al. 2012). 
Once the first stars are formed, they shine and emit radi- 
ation. This can have relevant impacts (radiative feedback) 
(Efstathiou 1992) on the surrounding medium and on the 
following star formation history (e.g. Ciardi et al. 2000, 2001; 
Whalen et al. 2004; Iliev et al. 2005; Alvarez et al. 2006; Susa 
& Umemura 2006; lUev et al. 2006, 2009; Whalen et al. 2010; 
Mellema et al. 2006; Susa et al. 2009; Gnedin et al. 2009; 
Petkova & Springel 2011; Paardckoopcr et al. 2011; Wolcott- 
Green et al. 2011). Therefore, appropriate radiative transfer 
(RT) calculations coupled with hydrodynamics and chem- 
istry must be performed in order to consistently explore the 
repercussions on molecule destruction or enhancement, and 
hence on structure formation in the early Universe. There 
have been studies on how the radiation feedback affects 
larger halos such as clusters and galaxies (e.g. Shapiro et al. 
2004; Yoshida et al. 2007; Croft & Altay 2008), and par- 
ticular interest has arisen from popHI stars (Wise & Abel 
2008; Hasegawa et al. 2009; Smith et al. 2009: Whalen et al. 
2010), whose large luminosities (which suggest them as pro- 



genitors of early gamma-ray bursts, e.g. Campisi et al. 2011; 
de Souza et al. 2011; Campisi et al. 2011) not only ionize, but 
also expel gas from the pristine mini-halos they sit in. For 
example, Lyman- Werner photons (11.2 — 13.6 eV) at early 
times could eradicate H2 from halos, delaying or completely 
impeding the collapse of molecular gas (e.g. Susa 2007; Wise 
& Abel 2007; .Johnson et al. 2007; Ahn et al. 2009; Trenti 
& StiavelU 2009; Whalen et al. 2010). Non the less, several 
issues in this respect are still unsolved, like the role of ra- 
diative feedback on popIII star formation, its efficiency in 
halting or boosting gas cooling in primordial environment, 
its interplay with mechanical and chemical feedback, its ef- 
fects on the dynamical and thermodynamical state of the 
cosmic gas, or its connections to the set-up of a turbulent 
medium dominated by hydro-instabilities (e.g. Maio et al. 
2011). 

First semi-analytical works (Haiman et al. 1997a,b; Haiman 
1999) expected molecular hydrogen to be universally de- 
stroyed by UV stellar radiation in the LW band. 
Subsequent, one-dimensional, numerical studies of small- or 
intermediate-size boxes (< 1 Mpc a side), including RT cou- 
pled with hydrodynamics, demonstrated that these predic- 
tions were overestimating molecular destruction, and, de- 
spite the stellar UV radiation, H2 could be rapidly reformed 
in the shock compressed gas of the ionization fronts (I- 
fronts) of HH regions (Ricotti et al. 2001, 2002a,b, showed 
that via a Softened Lagrangian Hydrodynamics Particle- 
Particlc-Particle-Mesh, SLH-P''*M, implementation). This 
was the first suggestion for "positive" feedback exerted by 
I-fronts on primordial structure formation. 
Similar calculations applied to individual haloes nearby pri- 
mordial stellar sources (e.g. Shapiro et al. 2004; Iliev et al. 
2005; Hasegawa et al. 2009) showed that UV radiation from 
popIII stars in crowded star forming regions could photoe- 
vaporate small (~ 10^ M©) haloes in > 100 Myr. However, 
more recent studies (e.g. Susa & Umemura 2006; Susa et al. 
2009; Hasegawa et al. 2009; Whalen et al. 2008; Whalen & 
Norman 2008; Whalen et al. 2010) found that halo photoc- 
vaporation due to first stars strongly depends on the features 
of the stellar sources and on the hydrodynamical properties 
of the collapsing cloud. Thus, radiative feedback was found 
to be able to trigger gas cooling by catalyzing H2 formation. 
In order to overcome numerical issues in larger cosmological 
simulations (> IMpc a side), further studies adopted sim- 
ple analytical prescriptions for a mean, uniform, UV back- 
ground, assumed to be established after the early onset of 
star formation (e.g. Machacek et al. 2001, 2003; Mesinger 
et al. 2006, 2009, via a grid Adaptive Mesh Refinement, 
AMR, implementation). In the same spirit, detailed, AMR, 
popIII, star formation calculations in uniform LW back- 
grounds (O'Shea & Norman 2008a; Wise & Abel 2008; 
Shang et al. 2010) confirmed that the LW background is 
much less destructive to new popIII stars than originally 
supposed, and, even by assuming extremely high background 
values of J ^ lO"^'^ ergs~^ cm~^ Hz"'^ sr~^ (= J21) star for- 
mation is just delayed, not shut down. 
This essential review shortly shows how debated the issues 
related to radiative feedback axe and highlights the numer- 
ous unknowns when dealing with RT and molecular chem- 
istry. In the present work, we aim at contributing to the sci- 
entific discussions and works already existing in literature 
by implementing three-dimensional RT schemes fully cou- 
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pled with hydrodynamics, and non-equihbrium chemistry. 
At the end of the paper, together with chemistry and ra- 
diative transfer, we will also consider additional feedback 
mechanisms which are of relevant interest for cosmological 
structure formation, namely, SN feedback and wind food- 
back. This represents a step further with respect to previous 
implementations, since they usually focus on individual feed- 
back processes. We perform the implementation within the 
widely used and well-tested numerical N-body/SPH code 
GADGET (Springel 2005), in its most recent and updated 
version. In this way, we will be able to self-consistently study 
the effects of RT on chemical abundances, and to pinpoint 
the basic consequences for the cosmological evolution of the 
structures in the Universe. 

In Sect. 2 we describe the implementations of radiation 
(Sect. 2.1) and cosmic chemistry evolution (Sect. 2.2). In 
Sect. 3, we test our implementation by performing analy- 
ses of the Stromgren sphere problem (Sect. 3.1) in a static 
(Sect. 3.1.1) and dynamic (Sect. 3.1.2) density field, and 
chemical abundance evolution (Sect. 3.2). We then apply 
our method to cosmological structure formation simulations 
(Sect. 4). We summarize, discuss and conclude in Sect. 5. 



2 IMPLEMENTATIONS 

We use the parallel tree N-body/SPH (smoothed particle 
hydrodynamics) code GADGETS, an extended version of the 
publicly available code GADGET2 (Springel 2005), and mod- 
ify it in order to couple chemistry evolution and RT. In the 
following paragraphs, we give the details about the imple- 
mentations of the RT (Sect. 2.1) and chemistry (Sect. 2.2) 
parts, and, in the next section, we will show results from our 
test runs. 



2.1 Radiative transfer 

To follow the propagation of radiation we use the implemen- 
tation of RT in GADGETS (Petkova & Springel 2009) and 
expand it to a multi-frequency scheme. The original imple- 
mentation is based on a moment method, where the closure 
relation is the optically thin variable Eddington tensor, sug- 
gested by Gncdin & Abel (2001). To follow the transport of 
radiation, we solve the equation of anisotropic diffusion for 
the photon number density per frequency n-y{i'): 



dn.,(u) _ d ( 1 dn^iv)h'' 
dt dxj I /t(i') dxi 



where t is time, Xi and the coordinate components, 

c is the speed of light, k(z/) is the absorption coefficient, ft'^ 
are the components of the Eddington tensor, s^{iy) is the 
source function, and the Einstein summation convention is 
adopted for all exponents i and j. 

The Eddington tensor is obtained by summing up the 
contributions from the sources of ionizing photons and its 
components are given by 



p^3 



Tr(P'. 



where 

P'^(x) 



'dVp.(x') ("-f:^;,-/'^- 



(2) 



(3) 




hp!/ [eV] 

Figure 1. Intensity of a black-body spectrum as a function of 
photon energy in eV for two different effective temperatures - 
S X 10'* K and 10"'' K. The vertical line marks the end of the Lyman- 
Werner band at IS. 6 eV. 



is the radiation pressure tensor and p* is the stellar density. 

The tensor is computed via a tree and effectively removes 
the dependence of the scheme on the number of ionizing 
sources - an advantage in cosmological and galaxy formation 
simulations. 

The source term is treated in a Strang-split fashion, 
where photons are first "injected" into the medium sur- 
rounding the sources, and are then "diffused" via equa- 
tion (1). Solving the equation for all particles in the sim- 
ulation reduces to a linear system of equations, which we 
solve implicitly by using a Conjugate-Gradient scheme, that 
ensures robustness and stability even for large timesteps. 

For our multi-frequency extension we need to transform 
the photon number density to ionizing intensity and vice 
versa. Since the photoheating and photoionization rates in 
equations (7) (discussed more in detail in Sect. 2.2) are ob- 
tained by integrating the intensity over frequency, we use a 
single photon number density in each frequency bin. 

The photon number density per frequency is derived 
from the ionizing intensity, /(!>"), as 



1 inlju) 

n-tiv) = r 

c hvV 



(4) 



where 4tt is the full solid angle and hp is the Planck constant. 
For any particle and at any timostop, the RT equation (1) 
is then solved for each frequency bin. 

The intensity of the radiative source, /(f), depends on 
the particular problems treated. For stellar sources, a com- 
mon, simple and suitable approximation (which we will also 
adopt in the following) is a black-body spectrum (see Fig. 1), 
with effective temperature dependent on the assumed stel- 
lar population. This treatment allows us to take into account 
contributions from a wide frequency range, even below the 
H-ionization energy of 13.6 eV. Indeed, as evident from the 
plot in Fig. 1 the low-energy tail of the black-body spectrum 
contributes to the Lyman- Werner band and that will defi- 
nitely affect the molecular evolution of the gas (see later) . 

We note that dealing with radiation below ^ 13.6 eV is a 
very debated and complicated problem (e.g. Thoul & Wein- 
berg 1996; Haiman ct al. 1997a,c; Omukai & Nishi 1999; 
Machacek et al. 2001; Kitayama et al. 2001; Ricotti et al. 
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Table 1. Ionization energies in [eV] 



H H- He He+ D H2 HD HeH+ 
13.6 0.7 24.6 54.4 14.9 15.4 15.4 44.5 



2002a,b; Mackey et al. 2003; Shapiro et al. 2004; Dijkstra 
et al. 2004; Susa & Umemura 2004; Stacy et al. 2011), be- 
cause of the many lines involved (76) in the LW band. Pull, 
detailed modeling of RT in such regimes is beyond the aims 
of this work, and, at some levels, it might be superfluous (e.g. 
Ricotti et al. 2001), since the ionizing flux at high redshift 
(z) is dominated by radiation from neighboring haloes (e.g. 
Ciardi et al. 2000), rather than from the soft-UV background 
in the LW band. Thus, H2 photodissociation is simply ac- 
counted for by using (e.g. Abel ct al. 1997; Ricotti ot al. 
2001; Yoshida et al. 2003; Ahn & Shapiro 2007; Maio ct al. 
2007; Whalen & Norman 2008) the radiative rate obtained 
from integration of the source intensity over the LW range, 
11.2 - 13.6 eV - sec cq. (7), in the next Sect. 2.2 - shielded 
according to Draine & Bertoldi (1996). 

2.2 Chemistry 

To couple radiation with chemistry, we include non- 
equilibrium reactions for H, He and molecule evolution 
(whose ionization and dissociation energies are quoted in 
Tab. 1), by following a chemical network (see Tab. 2) of sev- 
eral species: e", H, H+, H", He, He+, He++, H2, H+, D, 
D"*", HD, HeH"*". Besides some updates in the rates and in 
the reaction network, the implementation used is the same 
as the one in Maio et al. (2007) and Maio et al. (2010) (and 
based on Abel et al. 1997; GaUi & Palla 1998; Yoshida ct al. 
2003). Since the main coolants at early times are H-derived 
molecules, H2 (e.g. Saslaw & Zipoy 1967) and HD (e.g. Lepp 
& Slmll 1984), the inclusion of a large network is crucial to 
correctly resolve the hydrodynamics and the fragmentation 
processes of high-redshift g cLS, clS well demonstrated by e.g. 
Abel et al. (1997), Yoshida et al. (2003, 2006), Maio et al. 
(2006, 2007, 2009, 2010, 2011), Maio et al. (2011); Maio & 
lannuzzi (2011). 

To take into account chemical evolution, at each 
timestep and for each species i, the time variation of its 
number density n, is computed, for collisional and photoion- 
ization/photodissociation events, via 

^ ^ ^ ^ ^ kpq^iTlpTlq ^ ^ huTllTii hyiTli^ (5) 
P Q I 



where kpq^i is the rate of creation of the species i from species 
p and q, ku is the destruction rate of the species i from col- 
lisions with species /, and fc-yi is the photoionization or pho- 
todissociation rate of species i due to radiation (see Tab. 1). 
The collisional rates are given by 



-I' 



i{u)f{u)d u, 



(6) 



(and the analogous for ku), with u relative velocity of par- 
ticles p and q, cjpq^i{u) interaction cross-section, and f{u) 
MaxwcUian velocity distribution function. More precisely, 
we do not compute the integrals on the fly, but instead we 



interpolate pre-computcd tables (whose references are listed 
in Tab. 2) in order to speed-up the code. These rates are 
temperature dependent and are expressed in units of vol- 
ume per time, i.e. [cm^s"'^] in the cgs system. As in the 
left-hand side, the first and second term in the right-hand 
side of equation (5) have dimensions of a number density 
per unit time, [cm^'"' s^^] in the cgs system. 
Consistently with Sect. 2.1, when the species i interacts with 
radiation (7) - see Table 2 - and gets photo-ionized (like for 
H, D, H^, He, He+) or photo-dissociated (like for H2, Hj, 
HD, HeH''"), the corresponding radiative rate k^i can be 
written as 



(7) 



where the An stands for isotropic radiation, is the 

source intensity as a function of frequency v, (T.yi{i') is 
the cross-section for the given process, hp is the Planck 
constant, c the speed of light, and n^f{i') the photon number 
density per frequency. In the final equality we made use 
of eq. (4), to formally get the rate expression similar to 
eq. (6). The radiative rates are probabilities per unit time, 
and are given in [s~^] in the cgs system. So, the number 
of radiative interactions per unit time and volume between 
photons and Ui particles The latter quantity is 

added on the right-hand side of equation (5) when photon 
interactions are taken in consideration, and consistently 
with the other terms in the equation, this is also given in 
units of number density per time, i.e. [cm^"^ s^^] in the cgs 
system. 

When talking about radiative interactions, one has to 
consider that, while ionization energies (see Tab. 1) are 

uniquely defined, molecular dissociation energies might 
depend on the particular radiative process considered, and 
thus all the various channels must be taken into account. 
For example (see Tab. 2), LW radiation above the energy 
threshold of 11.2 eV can dissociate H2 in 2H, but harder 
photons with energies larger than 15.42 eV would simply 



ionize the residual H2 into Hj -|- e 



Also for HJ there are 



two possible branches: one above an energy threshold of 
2.65 eV and below 21 eV (HJ -|- 7 H+ -|- H), and a second 
one between 30 eV and 70 eV (H+ -I- 7 ^ 2H+ -|- e"). For 
the HeH"*" radiative interaction considered in our network 
the energy threshold is about 1.7 eV. 

The set of differential equation (5) is integrated via sim- 
ple linearization, so, given the timestep At, at each time t 
the temporal variation of the number fraction of species i 
can be written as 



t-HAt t 

[H [H_ _ Qt+At _ jjt+At t+At 

At ' i i ' 



(8) 



where we have introduced the creation coefficient for the 

species i 



— ^ ] y ] kpq^iUpUq, 



(9) 



p 9 

-3„-ll 



in [cm s ], and the destruction coefficient 

Di = ^ kiifii + k-yi, (10) 

I 

in [s~^]. The contribution from photoionization or photodis- 
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sociation is accounted for by adding, in equation (10), the 
k^i rates. The number density is updated from equation (8): 



Table 2. Reaction network 



t+At 



1 + D. 



t+At 



At 



(11) 



We apply this treatment to all the chemical species included, 
with the coefficients for each reaction in the network quoted 
in Table 2. Gas coohng or heating is computed from H and 
He collisional excitations (Black 1981; Cen 1992), ionizations 
(Abel et al. 1997), recombinations (Hui & Gnedin 1997b), 
H2 and emissions (Galli & Palla 1998), HD emissions 
(Lipovka et al. 2005), Compton effect, and Bremsstrahlung 
(Black 1981). 

The timestepping is limited by the cooling time, 
E 

1 ' 



tcool — 



(12) 



and by the electron recombination time. 



te 



(13) 



where E and rie are the energy and the electron number 

fraction of each particle, and E and he the corresponding 
time variations. For accuracy reasons in the abundance de- 
terminations, the chemical subcycles are done over 1/10 of 
min(tcooi,te) (as described by e.g. Anninos et al. 1997; Abel 
et al. 1997; Yoshida et al. 2003; Maio et al. 2007). The addi- 
tional constrain given by te is useful mostly when rie changes 
very steeply, like behind the shock fronts, or at ~ 10* K, be- 
low which hydrogen recombines very efficiently and above 
which hydrogen gets ionized very fast. Further details can 
be found in Maio et al. (2007) and references therein. 



3 TEST SIMULATIONS 

In order to test the implementation we perform numeri- 
cal simulations under different conditions. First, we numer- 
ically solve an expanding ionized sphere problem (Sect. 3.1) 
by fully including both the RT and chemistry treatments. 
Then wc show the cosmic evolution of the different chemical 
species, coupled with the radiative gas emissions (Sect. 3.2). 
Finally, we will perform cosmological simulations of early 
structure formation (Sect. 4) to check the effects of radia- 
tive feedback on gas cooling and collapse. 

For all simulations wc use a set of 284 frequencies, cov- 
ering the range from 0.7 eV to 100 eV. The density of fre- 
quency bins around the peaks of the photoionization and 
photodissociation cross-sections is increased - i.e. there are 
more frequency bins in the spectral regions of interest. In 
this way wc can ensure that photons arc both traced and 
absorbed properly and no spectrum-averaging mistakes are 
made. 



3.1 Ionized sphere expansion 

3.1.1 Ionized sphere expansion in a static density field 

The expansion of an ionization front in a static, homoge- 
neous and isothermal gas is the only problem in radiation 
hydrodynamics that has a known analytical solution and 
is therefore the most widely used test for RT codes (e.g. 
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+ H 2H + e~ 


A97 / Y06 / M07 


+ H+ — > 2H 


P71 / GP98 / Y06 / M07 


+ H+ — > H2+ + c~ 


SK87 / Y06 / M07 


H2+ + ^ 2H 


GP98 / Y06 / M07 


TU +- 1 TU — ( TU 1 lU 

tl2 -|- tl — ^ tl -|- n2 


A AT / ID AO / "VAC / A /TAT 

Ay7 / GFyo / YUd / MU7 


H2 + 7 H+ + e- 


A97 / YOd / M07 


-n.2 "r 7 — r z n 


PCd( 1 xlUl / I Uo / iViUi 


D -(- H2 HD -1- H 


WS02 / M07 


D+ -1- H2 ^- HD + H+ 


WS02 / M07 


HD -1- H ^. D -1- H2 


SLP98 / M07 


HD -1- H+ ^. D+ -1- H2 


SLP98 / M07 


H+ + D -i. H + D+ 


S02 / M07 


H -1- D+ -i. H+ + D 


S02 / M07 


D+ -1- e- D -1- 7 


GP98 


D -1- 7 ^ D+ -(- e- 


GP98 


He -f H+ -> HeH+ + 7 


RD82/ GP98 / M07 


HeH+ + H ^ He + H+ 


KAH79 / GP98 / M07 


HeH+ -(- 7 -> He -f H+ 


RD82 / GP98 / M07 



Notes: 7 stands for photons; P71 = Peterson et al. (1971); 
KAH79 = Karpas et al. (1979); RD82 = Roberge & Dalgarno 
(1982); SK87 = Shapiro & Kang (1987b); A97 = Abel et al. 

GP98 = Galli & Palla (1998); SLP98 = Standi et al. 
ST99 = Stibbe & Tennyson (1999); ROl = Ricotti 
et al. (2001); WS02 = Wang & Stancil (2002); S02 = Savin 

(2002) ; GB03 = Glover & Brand (2003); Y03 = Yoshida et al. 

(2003) ; S04 = Savin et al. (2004); Y06 = Yoshida et al. (2006); 
M07 = Maio et al. (2007). 



(1997) 
(1998) 



Iliev et al. 2006, 2009). For such a set-up, the ionized bub- 
ble around the ionizing source reaches a final steady radius, 
called the Stromgrcn radius, where absorptions and recom- 
binations are balanced along the line of sight. For an H-only 
gas, the Stromgren radius is analytically given by 



= (- 



3iV^ 



Y 47raBng 



1/3 



(14) 



with N~f - the luminosity of the source in photons per sec- 
ond; ttB - the case-B recombination coefficient; and nn - 
the hydrogen number density. The case-B recombination co- 
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efficient assumes the so called 'on-the-spot' approximation, 
where photons from recombinations to lower energy levels 
are immediately absorbed in the vicinity of their emission 
(e.g. Spitzer 1978). If we approximate the ionization front (I- 
front) as infinitely thin, i.e. it features a discontinuity in the 
ionization fraction, the temporal expansion of the Stromgren 
radius can be solved analytically in closed form, with the I- 
front radius ri given by 



rs[l- 
where 



■ exp{-t/trcc)] 



1/3 



^rcc — 



(15) 



(16) 



is the recombination time. 

In our first test we perform an ionized sphere expan- 
sion, but we allow the temperature of the gas to vary in 
order to test the coupling between the RT and the full non- 
equilibrium chemistry treatment. As a reference, we com- 
pare to the analytical case with constant temperature. We 
follow the expansion of an ionized sphere around a source 
that emits iV.^ = 5 x 10** photons s~^. The shape of the 
source spectrum corresponds to a 3 x 10** K black body. 
The surrounding gas density is p = 1.7 x 10^^'^gcm^'^ 
(~ 10"^ cm"^) and is sampled by 16^, 32^, and 64^ gas par- 
ticles^ . In the 32'^ case also the shielding of Draine & Bertoldi 
(1996) has been adopted, with the values cited in their pa- 
per. The initial temperature of the gas is set to T = 10^ K 
and is subject to photoheating and radiative cooling. At a 
temperature of 10* K, the case-B recombination coefficient is 
qb = 2.59x10"^^ cm^s"^ (e.g. Iliev et al. 2009). Given these 
parameters, the recombination time is t^cc ~ 125.127 Myr, 
and the expected Stromgren radius in the isothermal case 
(assuming T = 10* K) is rs — 5.4 kpc. In Fig. 2, we show 
the evolution of the radial position of the I-front with time 
for the different resolutions. As a proxy for the position of 
the front we take the radius where the neutral and ionized 
hydrogen fractions are equal (see also Fig. 3). 
All resolutions agree very well with each other. There is 
no difference in the case with shielding since the simula- 
tion never reaches the densities required to produce some 
effect, as discussed in the introduction. Our results agree 
within 10% with the analytical ones from equation (15). 
In particular, the simple analytical solution is systemati- 
cally larger than the full-simulation trend, which is expected. 
This can be explained by the missing cooling contributions 
in the analytical calculations from, e.g.. He, II2, Hj, HD 
that lower temperatures, enhance recombination, and make 
the Stromgren radius decrease (as visible in the simulated 
case). In fact, equation (15) is computed by assuming con- 
stant temperature for hydrogen-only gas (see also Petkova 
& Springel 2009), while, in the numerical calculations the 
full chemistry treatment of Table 2, including cooling and 
heating, is considered. Pawlik & Schaye (2011) find similar 
results in their one-dimensional ionized sphere simulations^ . 



^ For more resolution studies, see Petkova & Springel (2009). 
There, it is shown that numerical convergence is reached already 
with 8'^ particles. 

^ Note that Pawlik & Schaye (2011) adopt a black body spec- 
trum with effective temperature of 10^ K and therefore produce 
different temperature and He radial profiles. 




200 300 
t [Myr] 



Figure 2. Evolution of the radial position of the I-front. The 
red line shows the analytic solution for an isothermal hydrogen 
only sphere. The black lines show our results from the simulations, 
assuming the radius is at the position where the amount of neutral 
and ionized hydrogen is equal. The thick dashed line shows the 
results from the simulation with 64^ particles, the thinner dash- 
dotted line — 16^ particles, and the thin solid line - 32^ particles, 
where also H2 shielding has been adopted. All lines agree very well 
with each other, where the lowest resolution exhibits more scatter. 
The results from the simulation with shielding show that it can 
be discarded in set-ups like this, where no high particle number 
densities are reached. As shown in other studies (see text), the 
radius of the I-front stays always below the analytical solution 
given by eqn. 15. 




Figure 3. Top panel: Chemical abundance fractions radial pro- 
files at 500 Myr after the source has been switched on. Bottom 
panel: Temperature radial profiles at 500 Myr after the source 
has been switched on. The temperature inside the ionized region 
reaches ~ 10* K and extends beyond 6 kpc since harder photons, 
unabsorbed by the gas, heat the medium ahead of the ionization 
front. All simulations with different resolution ad shielding agree 
very well with each other at all radii. 
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In Fig. 3, we show the radial profile of the temperature 
of the gas at 500 Myr after the source has been switched on. 
The temperature inside the ionized region reaches ~ 10'' K, 
consistently with photoheating from a stellar-type source, 
and extends beyond 5 kpc. Even further the temperature 
begins to drop. Harder photons (with energies > 60 eV), 
that do not ionize the elements effectively, heat the medium 
ahead of the ionization front. If the source had a harder spec- 
trum, e.g. 10^ K black body, than the gas would be heated 
even at larger radii. We stress that at a distance of > 5 kpc 
(namely, around the Stromgren radius) temperatures steeply 
drop from ~ 10* K down to ~ 10^ K and recombination pro- 
cesses take place (see next). 

Correspondingly to the temperature profile, in Fig. 3, 
we display the radial profiles of the different chemical 
abundances at 500 Myr after the radiative source has been 
switched on. We assumed initial cosmic abundances^, which 
means that hydrogen species account for ~ 93% of the to- 
tal number densities and helium species for ~ 7%. As ex- 
pected, ionized fractions usually have larger values clo.sor to 
the sources (within a few kpc), while neutral or molecular 
fractions increase at larger distances (above 3 — 6 kpc) . In 
the following we discuss these trends, by referring to the re- 
action network of Table 2, in a more precise and detailed 
way. 

• Atomic hydrogen species: due to the strong radi- 
ation intensities near the central source, hydrogen is kept 
completely ionized (i.e. the total hydrogen fraction, ~ 0.93, 
is in the ionized state, H"^) within r < 5 kpc by the dominant 
photoionization process 



H-h7 ^ H+ -^e" 



(17) 



Only at larger radii, where the radiation intensity decreases, 
the recombination process 



H+ -h e" H -h 7 



(18) 



takes over and makes the H"*" fraction drop down of sev- 
eral orders of magnitude, with the consequent increase of H 
fraction up to ~ 0.93 (namely, the total hydrogen fraction 
is now in the neutral state). 

is a very important species because it represents one 
of the main channels by which molecular hydrogen can be 
formed (see discussion later). It shows fractional values of 
~ 10~® — 10"^° until r < 6 kpc and at larger distances drops 
of roughly 3 — 4 orders of magnitude. The increment of H~ 
fraction up to ~ 10^® in correspondence of r ~ 4 — 5 kpc is 
due to the simultaneous hydrogen recombination - eq. (18) 
- which leads to an increase of the neutral hydrogen and 
enhances the H~ formation via 



H-l-e ^-H -1-7 



(19) 



reaction. In the innermost regions, due to the low ionization 
energy of only 0.755 eV, H~ photoionization 



H -1-7 ^- H-l-e 



(20) 



3 



They are set to: x^- ~ 4 x 10"*, xn = 0.926, x^+ ~ 4 x 10""', 

= 10-19, XHe = 0.07, XHe+ = lO'^^, a:He++ = 10"=^°, 
— in— 18 — in— 16 ~„ — in— B ~ 



Xh- =-1-0 , XHe = U.Ul-, Xjje+ = iO a:He++ = 

= 10-13, = 10-18, xno = 10-i«, xu = 10-B, Xd+ 
10-'^, a;„„„+ = 10-21. 



and coUisional destruction by the abundant H"*" and e 
species 



H" -h H+ ^- 2H 
H" -h H+ ^- Ha -h e" 
H" -h e" ^- H -I- 2e" 



(21) 
(22) 
(23) 



determine a lower fraction of ~ IQ-i^. At larger distances 
(r ~ 5 — 10 kpc), H is dominant, but, contrary to reaction 
(19), free e~ are lacking and the most effective reactions are 



H" -h H Ha -h e" 



(24) 



H" -h H ^ 2H -h e" (25) 

that lead to a decrease of fraction down to < IQ-i^. 

We note that a calculation of the exact analytical expres- 
sion for the the Stromgren radius and a comparison with 
our results is not meaningful, since the analytical study is 
based on the simplified case of hydrogen-only gas and does 
not take into account interaction with other species and the 
effects of photoheating. However, the radius of the ionized 
hydrogen reaches ~ 5 kpc, which is approximate to the ex- 
pectation value of the Stromgren radius for the isothermal 
case. 

• Atomic helium species: in the same way, due to the 

reactions 



He -I- 7 He+ -|- e 
He -I- e" -> He+ + 2e" 



(26) 
(27) 



neutral He is efficiently destroyed at r < 3 kpc, and the 
residual fraction is ^ 0.01. He"*" and He''"'' reach fractions of 
~ 0.06 and ~ 10^'', respectively, via 



He+ + c" 
He+ 7 ■ 



He++ + 2e" 



He^ 



+ e" 



He'^'^ +e~ ^ He'^ + 7. 



(28) 
(29) 
(30) 

Moreover, the additional He depletion by H''' collisions leads 
to the formation of HeH''' (as we will discuss later). At 
r > 3 kpc, the decreasing intensity of ionizing radiation - 
which plays the most relevant role in equations (26) and 
(29) - and the ongoing recombination processes which take 
away free electrons from the medium - needed e.g. in reac- 
tions (27), (28) - cannot sustain ionization any longer and 
the trends exhibit monotonic radial drops for both He^ and 
He''"'' . More exactly. He''' and He recombine according to 



He 



++ 



+ e 



Ue^ -I- e" He -I- 7 



(31) 
(32) 



respectively, so, the final states of He''" and Hc^''' arc strictly 
linked to each other with abundances of ~ 10"^ and 10-" 
at r ~ 8 kpc. Radiative interactions are weaker and weaker 
and they become practically negligible at large radii. 

• Atomic deuterium species: the behaviour of D and 
D+ are quite regular and similar to H and H'*', with D''' 
being dominant at r < 5 kpc and D dominant at r > 5 kpc 
of ~ 3 orders of magnitude. The abundances of deuterium 
and hydrogen species are bound by the balance reactions 



''HeH+ 



H'^ + D ^ H -h D'^ 
H -h D'^ -5. H'^ + D, 



(33) 
(34) 



© RAS, MNRAS 000, 1-16 



8 M. Petkova & U. Maio 



but the most efficient processes are still photoionization 
('close' to the source) 



D + 7^-D++e" 

and recombination ('far' from the source) 
D+ + e" -)-D + 7. 



(35) 



(36) 



Minor contributions of D and D+ arc also involved for molec- 
ular species (sec next) and can slightly affect If or pro- 
duction. 

• Molecular species: Finally, we discuss the trends of 
molecular species (H2, Hj, HD, HeH"*"). Their profiles are 
less regular and intuitive than atomic profiles, since more 
processes need to be addressed at the same time. For ex- 
ample, H2 and have fractional values of ~ fO~^^ and 
~ 10~^°, respectively, at r ~ Ikpc, then, they show an in- 
creasing trend and a peak of ~ 10"^° — 10~^ at r ~ 4— 6kpc. 
For larger r, they exhibit a drop off below ^ 10~^*, with H2 
overcoming at r > 7kpc. These behaviours are under- 
stood by considering that in the innermost regions the ra^ 
diative negative feedback on molecule formation decreases 
for increasing r: this means that at larger radii radiation is 
not strong enough to dissociate molecules via 



Ha + 7 ^- 2H 



(37) 



and, thanks to the available H, and e~ (simultaneously 
present at r ~ 4 — 6kpc, where H recombination is still 
taking place), H2 formation can proceed through the H~ 
channel - see also reaction (19), 



H-fe" 



H 



-7, 



and through the channel, 
H + H+ ^ H+ 7 
+H^-H2-|-H+. 



(38) 
(39) 

(40) 
(41) 



This is the intrinsic reason why molecular hydrogen roughly 

follows the H~ profile discussed earlier. Obviously, residual 
photons will slightly boost the ionized fraction of (with 
respect to H2) via 



H2 + 7 ^- HJ + e- 



(42) 



until r ~ 7 kpc. At larger distances, photons are too weak to 
ionize the gas, so H2 takes over and drops dramatically 
of ~ 3 orders of magnitude within 1 kpc. Additionally, we 
note that the sharp decrement of molecular fractions at very 
large distances is basically due to the fact that the medium 
becomes almost completely neutral and the ionized fractions 
of e~ and H"*" are too low to boost H~ and HJ, and, hence, 
H2 formation. 

Similarly, the increase with radius of HD fraction is es- 
sentially caused by the weakening of the central radiation 
and of H^ and c^ fractions which are less and less effective 
in dissociating it at larger r. In particular, at r < 4 kpc, 
H, D and H2 are subdominant with respect to their ionized 
counterparts H"*", D"*" and Hj, thus, HD formation due to 



D + H2 HD + H 



(43) 



is strongly inhibited, while little contributions come from 
D+-I-H2 -^HD-fH+. (44) 



Destruction from reaction 

HD H+ -> D+ H2 (45) 

further lower HD abundances around ~ 10^ level. At 

r > 4 kpc, hydrogen and deuterium recombinations, to- 
gether with H2 formation via H~ channel, support HD 
formation through reactions (43) and (44), instead reac- 
tion (45) is no longer effective. As a consequence, the HD 
fractional abundance grows more than ~ 2 orders of magni- 
tude at r ~ 8 kpc. 

We note that a boost of HD production in shock-compressed 
gas was expected because of the simultaneous increment of 
H2 fractions, quite visible in Fig. 3, in correspondence of the 
Stromgreen radius. 

For what concerns the aforementioned HeH"*", this is effi- 
ciently produced near the source because there are a lot of 
free protons which can boost its abundance via 

He H+ ^- HeH+ -|- 7 (46) 

even in a more powerful way than photodissociation 

HeH+-h7^He-|-H+. (47) 

Only when protons arc lacking (i.e. at large r) HeH^ pro- 
duction is inhibited up to ~ 3 orders of magnitude and drops 
from ~ 10"" to ~ 10"". 

We stress that the residual relative ionized fractions 



far from the source (r 



8 kpc) 



10" 



no+/nu ~ 10 ^, riH-ZriH ~ 10 riHe+Z^He ^ 10 

nHe++MHe ~ 10~^, n„+/nH2 ~ 10~^. The absolute values 

2 

quoted in the previous discussion are dependent on the 
initial composition assumed for the gas. Although the qual- 
itative behaviour is not supposed to change much, larger 
or smaller values for the assumed fractions could result in 
more efficient or less efficient formation and destruction 
processes. However, the general trends are supposed to be 
quite independent from that. 

The test presented here is based on the initial conditions 
of Iliev et al. (2009), but different works are available in lit- 
erature. Results that are going in the same direction were 
in fact found by e.g. Ricotti et al. (2001); Ahn & Shapiro 
(2007); Whalen & Norman (2008), who developed, indepen- 
dently, different H, He, and H2 chemistry networks. Their 
simulation setups were very different, though, and the pa- 
rameters for the central source, as well. 
Ricotti et al. (2001) (their Fig. 3) performed the first one- 
dimensional calculations of a star shining in a primordial 
mini-halo. A popIII-like source with a power-law spectrum 
and iV-y ~ 1.2 X 10"'^ at 2 ~ 19 was located in a static, 
uniform medium at the mean density n{z ~ 19) ~ 0.1 cm~^. 
They studied the behaviour of the molecular species after 
~ 100 Myr from the explosion of the central star, when the 
Stromgreen radius had reached ~ 4 kpc, and a correspond- 
ing peak of H2 with ~ 5 x 10^*, a few kpc wide, had 
formed. This was equivalent to a boost of about 3 orders of 
magnitude with respect to their initial ~ 10~^ value as- 
sumed. Also in our previous discussion, despite the different 
initial values, we have found an increase in the H2 fraction 
of ^ 3 orders of magnitude. 

Similarly, Ahn & Shapiro (2007) (their Fig. 8) and Whalen 
& Norman (2008) (their Fig. 1) performed one- and three- 
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dimensional calculations, respectively, considering a 120 M© 
popIII-like source with a blackbody spectrum having effec- 
tive temperature of Teff = lO'' K and N., ~ 1.5 x 10^" s"\ 
They assumed the star to be in a primordial halo with a 
truncated isothermal sphere density profile, whose central 
matter densities were reaching up>4xl0"^^gcm-^ They 
focused on the innermost core of the star forming regions, af- 
ter roughly 0.6 times the lifetime of the massive popIII stars 
(a few Myr), when the Stromgreen radius was still around 
< 30 pc and a corresponding peak of H2 ~ 10 pc wide had 
arisen. They also find an increase to XH2 ~ 10~^, a few or- 
ders of magnitude larger than the assumed initial value of 
XHa ~ 10"''. 

The first relevant difference with the study presented here 
is in the assumed A'^^ (our value is 5 X 10*** s"^): Ricotti 
et al. (2001) used an A^'^, 2.2 times larger, while Ahn & 
Shapiro (2007) and Whalen & Norman (2008) used a value 
300 times larger Consequently, the number of photons avail- 
able to ionize atoms or dissociate molecules is much bigger 
in their cases, than in ours. This explains why they find H2 
fractions more significantly destroyed in the inner regions. 
Then, molecule dissociation or creation is also strictly linked 
to the typical conditions of the medium. In the aforemen- 
tioned works, this is denser than in ours (of ~ 2 — 4 orders 
of magnitude), and have strong implications when comput- 
ing abundances, since higher densities allow easier coUisional 
dissociation in the inner, ionized regions. When tempera- 
tures fall below ~ 10** K, H2 creation becomes more efficient 
and larger free-electron number densities can lead to a no- 
ticeable increase of molecular fractions* . This is actually the 
reason why the profiles of our low-density gas in Fig. 3 show 
H2 increases up to only ~ 10"^", while very overdense gas 
can reach H2 fractions of ^ 10~* in front of the propagating 
I- front. 

We conclude by noticing that, despite the huge diversities in 
the parameters and in the H2 peak values, the different be- 
haviours of the profiles are in qualitative agreement, showing 
an increment of the molecular fractions in correspondence of 
the I- front, and orders-of-magnitude declines farther away. 
Together with the H2 raise (within a shell of ~ 2kpc) the 
consequent HD increment is obviously expected to be more 
prominent in high-density environments. 
None the less, the final effects on triggering star formation 
are not completely clear, since gas runaway collapse can be 
significantly enhanced only when densities are larger than 
~ 1 — 10 cm~'^, and molecular fractions increase up to > 10~^ 
(see further discussions in Sect. 4). 




100 200 300 400 500 

t [Myr] 



Figure 4. Evolution of the radial position of the I-front. The 

dashed line shows our results from the simulation, assuming the 
radius is at the position where the amount of neutral and ion- 
ized hydrogen is equal. The solid line is the result from eq. (15), 
computed by assuming isothermal gas at 10* K. The extension of 
the solid line is given by eq. (48). The results from the simulation 
agree very well with the analytical results, also at radii beyond 
the Stromgren radius rg. 



t = 200 Myr 






3.1.2 Ionized sphere expansion in a dynamic density field 

In this section we repeat the test from the previous one, but 
we allow the gas particles to move due to pressure forces, 
but not gravity. The density field is not static anymore. The 
resolution we have chosen here is 32'^ gas particles. 

The position of the I-front in this stage is given by 
(Spitzer 1978) 



Figure 5. Top panel: Chemical abundance fractions radial pro- 
files at 500 Myr after the source has been switched on. Bottom 
panel: Temperature, density, and Mach number radial profiles at 
200 Myr after the source has been switched on. The tempera- 
ture inside the ionized region reaches ~ 10* K and extends be- 
yond 6kpc since harder photons, unabsorbed by the gas, heat the 
medium ahead of the ionization front. A shock, a H2 shell, and 
an HD boost have developed ahead of the I-front. 



* For example, electron fractions of ^ 10~* determine typical 
IGM fractions for H2 of ~ 10~®, instead, in gas with overdensities 
of ~ 10*, H2 fractions can be boosted up to ~ 10~^ (Ahn &; 
Shapiro 2007). 
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Figure 6. Number fraction evolution as a function of redshift for 
tlie different chemical species, followed in our implementation. 
The background cosmology is a standard flat ACDM model with 
geometrical parameters: f2o,tot = 1-0, f!o,A = 0.7, r2o,m = 0.3, 
!^o,b = 0.04. 

where Cs is the sound speed of the ionized gas and rs is the 
Stromgren radius given by equation (15). 

In Fig. 4 we show the evolution of the position of the I- 
front and its velocity with time. We compare with analytical 
results from equations (15) and (48). We note that the I- 
front follows the analytical prediction for static gas in the 
beginning of the expansion. After approximately 200 Myr at 
approximately 5 kpc from the source, the I-front continues to 
move outwards, rather than start to decelerate, and follows 
the dynamic gas solution given above. 

In Fig. 5 we present the radial profiles of the differ- 
ent species, the temperature, and the density of the gas at 
t = 200 Myr after the source has been switched on. The 
temperature in the ionized sphere is ~ 10'* K and there is a 
shock and a density contraction ahead of the I-front, which 
moves at super-sonic speed at this evolution time. We also 
note that there is a HD boost and a H2 shell ahead of the 
I-front as well, inside the density contraction region. Similar 
results have been shown also by Ricotti et al. (2001). These 
conditions - increased density and increased H2 fraction are 
favorable for star formation. 

3.2 Cosmological abundance evolution 

In order to test how our implementation performs during 
cosmological evolution, we run the chemical network cou- 
pled with the RT network for the mean background density 
evolution and follow the changes in the different species as 
a function of the cosmic time. The radiative source is as- 
sumed to be the uniform CMB radiation with a black body 
spectrum with effective temperature of ~ 2.73 (l + z) K. The 
CMB radiation is self-consistently followed, according to the 
treatment outlined in Sect. 2. The initial fractions, x, for the 
different species^ are initialized, at redshift z ~ 10'^, accord- 



ingly to a neutral plasma at ^ 10^ K (see e.g. Galli & Palla 
1998). 

We present the results in Fig. 6, where the cosmic mean 
number fractions as a function of redshift are plotted for a 
standard flat ACDM cosmology with geometrical parame- 
ters fictot = 1.0, D.o,A = 0.7, D.o,m = 0.3, ^o.b = 0.04. The 
H and He number fractions are unaffected by the redshift 
evolution. 

At early times, some residual recombination process 
continue taking place, while the CMB temperature goes 
down, and make H""" evolution drop from a fraction of 
~ 10"^, at 2 ~ 1000, to the final < 10"" value. The evolu- 
tion of H~ clearly shows the effects of free electrons at high 
redshift that boost its formation (and H-derived molecule 
formation) from Xjj- ~ 10~^°, at z ~ 1000, to x^- ^ 10~^^, 
at 2 ~ 100. The following hydrogen recombination implies a 
decrease both in H""" and e~ fraction, with consequent drop 
of H~ (see Table 2) down to x^- ~ 10~^^ at low redshift. 
The two ionization states of He are constantly kept at very 
low values, close to the initial ones (their trends are not 
displayed in Fig. 6 for sake of clarity). 

As already mentioned, hydrogen molecule formation is 
initially enhanced due to the available residual e~ and the 
progressively diminishing effects of CMB radiation: this al- 
lows hydrogen to increasingly form H^, with a peak around 
z ~ 300, and H2, until z ~ 70. Hj formation is efficiently 
driven by H and the available H^ in primordial times. At 
z < 300, paucity of free protons (whose fraction in the mean- 
time has dropped of about one order of magnitude) make 
difficult HJ formation. However, the continuous increment 
of H~ enhances H2 at 2 ~ 100, and afterward its formation 
is mainly driven by the H~ channel, rather than the Hj 
one, until z ~ 70. The dominant formation path at different 
times is clearly recognizable when comparing the H2 trend 
with the HJ and H~ trends. The z ~ 300 peak of Hj cor- 
responds to the steep increase of XH2 at early times, while 
the 2 ~ 100 peak of H~ corresponds to the following boost 
at later times. For 2 < 70, xh2 stays roughly constant be- 
tween 10"*^ — 10~^ because further production is halted by 
the decrement of free protons and electrons, and radiative 
destruction cannot take place because of the low CMB fiux 
at low redshift. 

Deuterium is also affected by cosmological evolution 
and the weaker CMB intensity eflects at lower redshift. 
Thus, recombination of D and D""" make xb increase up to 
- lO""* and Xn+ dramatically drop down at 2 < 100. The 
simultaneous ongoing H2 formation at 2 > 70 also 'drags' 
HD fractions up to ~ 10~^ levels (HD is very sensitive to 
H2 abundances). 

The HeH^ molecule is often formed behind fast shocks 
(e.g. Neufeld & Dalgarno 1989) by He and H+, but, due 
to its low dissociation energy (of ~ 14873.6 cm~^ ~ 1.7eV; 
Bishop & Cheung 1979) it can be found and emitted only 
below ~ 10* K. Moreover, the presence of background radi- 
ation can dissociate it in its two components soon after the 
recombination epoch, at 2 < 10"^. Indeed, the plot in Fig. 6 
highlights how HeH""" is initially underabundant and then is 
gradually formed while the Universe expands, cools and the 



5 



They are set to: x^- ~ 4 X 10^4, ^.^ = 0.926, ^ 4 X lO"", =10 =10 xhd = 10 xd = 10 5, = 

fj_ = 10-19, XHc = 0.07, Xh,+ = 10-25, Xh,++ = 10-30, 10-^a;H^H+ = ^-^^ . 
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CMB radiation gets weaker. Over the cosmic time a;HGH+ 
slightly increases of one order of magnitude, from ~ 10~^* 
up to > 10~^^, and, due to H interactions, it sustains 
with subsequent H2 formation. 

The low radiation intensity of the CMB in not able to 
produce large changes in the abundances of the elements, 
but is never the less a good test on the performance of our 
implementation. Our results agree very well with previous 
cosmological abundance evolution studies as e.g. Abel et al. 
(1997); GaUi & Palla (1998); Maio et al. (2007). 

The presence of an additional cosmic UV radiation dur- 
ing reionization (e.g. Haardt & Madau 1996) at low redshift 
(a < 7), would heat the medium and change the ionization 
equilibria. As a consequence, one might expect more free 
electrons, larger and abundances and hence more 
H2 production, accompanied by increased D"*" and HD frac- 
tions, and dissociation of HeH"*". 



4 APPLICATION: COSMOLOGICAL 
STRUCTURE FORMATION 

As a final application, we present results from a cosmolog- 
ical structure formation simulations in the framework of 
the standard ACDM model, and we will also consider ad- 
ditional feedback mechanisms from star formation (as al- 
ready mentioned in Sec. 1). We use a periodical comoving 
boxsize of Lbox = 0.5/i~^Mpc on a side with 2 x 64^ gas 
and dark-matter particles (for a resulting spacial resolution 
of ~ 0.4kpc//i comoving), sampled at the initial redshift 
z = 100. The runs include gravity, hydrodynamics, wind 
feedback, non-equilibrium chemistry, and radiative transfer. 
Stars are taken to be sources of ionizing radiation. Since 
each star particle in the simulation represents a whole stel- 
lar population with a Salpeter distribution, about 12% of its 
mass is in high-mass stars (> IOM0) that are able to pro- 
duce UV photons. The frequency distribution is assumed to 
be a black-body spectrum with an effective temperature of 
3 X 10'' K, which corresponds to a luminosity of approxi- 
mately 8 X 10*^ photons s~^ per high-mass star in the stellar 
population of the star particle. We assume that gas particles 
are converted into stars once a critical density of ~ 10 cm~"^ 
is reached, and the gas temperature is below ~ 10*K to 
make sure that the gas is effectively cooling (see details in 
Maio et al. 2009). Star forming particles also experience SN- 
explosion feedback, which heats the gas above ~ 10'' K, and 
wind feedback, which expels gas with a typical velocity of 
~ 500km/s (see also Springel 2005, and references therein). 

A pictorial representation of the simulated box is given 
in Fig. 7, where we show mass-weighted temperature slices 
through the simulation volume at redshift z = 10.61 and 
z = 6.14 for the full run, including, in particular, both non- 
equilibrium chemistry and RT. 

The first sources are well visible close to the central part 
of the slice, and, while the structure growth proceeds, more 
sources are found in scattered places along the converging fil- 
aments. Obviously, in a wider perspective (say on ~ 100 Mpc 
scale), radiative sources would be much more uniformly dis- 
tributed. However, there are currently serious computational 
limitations for performing simulations with such large box 
sides and, simultaneously, with resolution good enough to re- 
solve chemical evolution and radiative transfer at the same 



time. In the maps, the filamentary cold structures led by 
early molecular gas are well visible at temperatures around 
hundreds Kelvin. In the densest regions, radiative effects 
from first stars heat the medium above ~ 10"* K already by 

redshift z > 10. More and more star formation episodes ap- 
pear at later stages and contribute to the cosmic reionization 
process down to redshift z ~ 6. 

Given the small size of our box, we can clearly focus on the 
infalling phases of the cold material in the intersection of pri- 
mordial filamentary structures, and on the subsequent SN 
explosions, which heat the gas and push material into the 
lower-density, void regions. This helps understanding the dif- 
ferent role played by the various feedback mechanisms, on 
the other side, the lack of a very large box also leads to 
an insufficient number of stellar sources to fully complete 
reionization by 2; ~ 6. 

In Fig. 8 we plot the phase diagrams (gas temperature 
vs. number density) for the simulations with and without ra- 
diative treatment. The colors refer to the II2 mass fraction. 
Molecular hydrogen is dissociated in the presence of ioniz- 
ing radiation (mostly in the Lyman- Werner band) , and this 
effect increases with decreasing redshift, due to the higher 
number of available photons. As a consequence, this low- 
ers the star formation process, since H2 provides the largest 
part of the sub-10'' K cooling. It is evident from the plots 
that without radiative feedback H2 can reach levels > 10~^ 
and boost early star formation. In presence of radiation from 
stellar sources, molecular fractions decrease of several orders 
of magnitude down to < 10^*. This effect is stronger for 
low-density gas, where molecules re-form more slowly. Such 
behaviour is clear from the both cases shown, at redshift 
z = 10.61 and 2 = 6.14. In the latter case, a wider tempera^ 
ture spread, due to thermal heating of the infalling gas at low 
density, appears as a consequence of the ongoing structure 
formation. Additionally, also SN-heated gas at 10^ — 10^ K 
is evident with extremely low molecular fractions (< 10~^®). 
The lack of very cold gas in the more realistic case with both 
chemistry and RT suggests that low-tomperaturo cooling by 
metals (e.g. Maio et al. 2007) could be an important mech- 
anism to sustain star formation after the first generation of 
stars (Maio et al. 2010)*^. To better underline the impacts of 
radiative effects on early chemistry, we also show in Fig. 9 
the redshift evolution of the mean and maximum number 
fractions of the two main molecules, H2 and HD, and of the 
basic species for the different channels of molecular forma- 
tion, H~ and HJ. 

As mentioned previously, after early star formation sets 
in, around redshift z ^ 12, the H2 fractions drops by several 
orders of magnitude (from average values of ~ 10~^ down 
to < 10"^°), as ionizing radiation starts to propagate in 
the simulated volume. This is not seen in the case without 
radiation, where H2 increases almost monotonically, with 
peak values of ~ 10"'^ at 2 ~ 6. 

The mean HD fraction stays roughly constant in the 
presence of ionizing radiation, but it is increasing in the 
case without radiation. When comparing peak values, one 
sees that, at 2 ~ 12, HD has fractions of ~ 10~^ in both 
cases, but in the RT case, it is significantly destroyed when 



^' One could also rely on reionized gas (Yoshida et al. 2007) to 
form subsequent baryonic structures. 
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Figure 7. Mass-weighted temperature slices through the box at redshift z = 10.61 (left) and z = 6.14 (right). The simulation includes 
feedback effects, full non-equilibrium chemistry and RT (see text). 





Figure 8. Phase diagrams - temperature versus gas comoving number density at redshift ^ ^ 10 (top) and redshift 2 ~ 6 (bottom) 
for a simulation without (left) and with (right) ionizing radiation. The color shows the mass fraction of H2. In the presence of ionizing 
radiation, H2 is depleted on all scales by many orders of magnitude. The effect is stronger at lower redshift due to the larger number of 
photons available. 



photons begin to propagate. Instead in the non-RT case its 
maximum values catch up with H2. 

The H~ fraction is crucial for H2 formation via the 
H~-channel: in presence of RT its mean values are lowered 
down to ~ 10~^^, so they cannot catalyze efficiently further 



molecule formation, while in the non-RT case, H would 
increase by two orders of magnitude. 

Similarly, catalyzes H2 formation via the HJ- 

channel, but photon propagation destroys molecules and in- 
hibits formation. As a comparison, in the non-RT run 
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Figure 9. Evolution of the mean and maximum H2, HD, H~ and fractions with redshift for the simulation with ionizing radiation 
(right) and without (left). When star formation sets in at redshift z ~ 11, the H2 fractions drop fast, where both the mean and the 
maximum values are affected. The mean HD fraction stays constant in the presence of ionizing radiation and is increasing with lower 
redshift in the absence of radiation. The maximum follows the same trend, with the difference that it drops down as the first photos 
begin to propagate. There is only a small change in the evolution of the H~ fraction. Finally, the fraction is slightly suppressed from 
the ionizing radiation. 



reaches mean values of ~ 10"^'', a couple of orders of 
magnitude larger than in the RT run. Maximum fractions 
for H~ and are similarly suppressed of about one order 
of magnitude. 

As a conclusion, the radiative feedback is responsible 
for molecule dissociation. Molecules are easily depleted not 
only from the external sources, but also from the central 
sources which have just been formed. This is justified by the 
fact that catastrophic molecular cooling sets in already at 
densities > lcm~^, when the material is optically thin, and 
thus there is no significant gas shielding preventing H2 and 
HD dissociation. In fact, in the simulations presented here, 
radiative feedback becomes effective at densities > lOcm"'^. 
In the innermost central regions (i.e. below a few hundreds 
of comoving parsec, not sampled by our simulations be- 
cause of resolution limits), densities could rise up to values 
greater than ~ 10* cm"'', become optically thick, and pro- 
duce significant shielding^. This should not affect the overall 
molecular destruction over ~ kpc scales by external sources, 
though. 

We stress that precise, quantitative assessments about the 
net effects of radiative feedback on the surrounding gas are 
difficult, because one should be able to prove RT on very 
different scales within N-body/hydro chemistry calculations. 
While destruction of molecules in low-density intergalactic 
medium (IGM) is intuitive and expected, the behaviour in 
the innermost cores of collapsed objects is quite debated 
and strongly dependent on the local density regime. In par- 
ticular, as pointed out by Susa (2007), the ignition timing 
of the source star is crucial in collapsing clouds: if the ig- 



^ If we consider that primordial haloes have radii of a few kpc, 
the fraction of volume interesting for such events is very small: for 
a protogalaxy having a radius of ~ kpc and developing a dense 
core within ~ 100 pc, the volume fraction of high-density, possibly 
shielded gas is < 10~^. 



nition takes place when the central density is larger than 
~ 10'' — 10*cm~^, the collapse cannot be halted by radia- 
tion. However, following in great detail the stellar-core for- 
mation and the ignition of nuclear reactions which lead the 
birth of the star is nowadays still prohibitive. 
As a consequence, cosmological simulations are inevitably 
affected by resolution issues, mostly at large densities, that 
limit the possibility of fully understanding how the various 
physical processes interplay in the very final stages of star 
formation. This also affects our capabilities of quantifying 
photon escape from the dense star forming regions. It is rea- 
sonable to think that in realistic conditions the inner cores 
will be optically thick and only after the photons will have 
heated up the cloud and thermalized with it they will be 
able to escape and dissociate molecules in the diffuse IGM. 
Definitive answers on this topics are still lacking, though. 
Furthermore, these many uncertainties even lead to differ- 
ences of several orders of magnitudes when comparing mass 
estimates of haloes affected by radiation (e.g. Ciardi et al. 
2000; Kitayama et al. 2001; Dijkstra et al. 2004). 



5 SUMMARY AND CONCLUSIONS 

Cosmic gas and star forming processes are fundamental key- 
words in modern Astrophysics. They are, however, very 
complicated to study because they lie in the highly non- 
linear regime. Thus, simple perturbative approaches do not 
give relevant information about gas collapse and baryonic- 
structure growth. To understand them better, it is crucial to 
consider all the involved physical and chemical mechanisms. 
In case of early structure formation, primordial chemistry 
has to be taken into account, since it is mainly via H-derived 
molecules that first objects can aggregate gas and form stars. 
Furthermore, it is also necessary to self-consistently include 
the RT of the photons produced by stellar sources. Indeed, 
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these travel into the cosmic medium and can ionize the sur- 
rounding gas or (hssociate formed molecules. As the impacts 
of RT on the chemical evolution of the early Universe are 
expected to be significant, it is important to couple such 
calculations in order to get a reliable picture. 
In this work, we have presented numerical methods of cou- 
pled chemistry treatment and multi-frequency RT, applied 
to N-body, hydrodynamical simulations. We used the SPH 
code GADGETS - an extended version of the publicly avail- 
able code GADGET2 (Springel 2005) - including the RT im- 
plementation of Petkova & Springel (2009), and the non- 
equilibrium chemistry implementation of Maio et al. (2007), 
following e", H, H+, H", He, He+, He++, H2, H+, D, D+, 
HD, HeH+. 

After a detailed description of the coupling between the 
RT treatment (Sect. 2.1) and the non-equilibrium chemi- 
cal treatment (Sect. 2.2), we have performed convergency 
tests (Sect. 3) and cosmological applications (Sect. 4) of our 
code. 

We have started (in Sect. 3.1) with the expansion of 

an ionized sphere around a stellar-type source emitting at 
a temperature of ^ 3 x 10 K in a uniform-density gas. We 
have traced the element and molecule evolution with time 
and compared to analytical results (in Fig. 2 and Fig. 3). 
We found that analytical analyses based on H-only gas at 
a constant temperature give a sufficient criterion to predict 
the evolution of the I-front. The different species reach dif- 
ferent ionization radii: we use the Stromgren radius defined 
for H-only gas at 10* K, as a reference. Deuterium species 
recombine at a radius comparable (within 10 per cent) to 
the Stromgren radius, while He gets completely neutral at 
^ 2/5 the Stromgren radius. H2 has an ionization radius 
that is typically larger than the Stromgren radius of about 
40 per cent (~ 7kpc vs ~ 5kpc). The compression of gas 
due to the propagating I-front enhances molecular species 
(H2 and HD) in the outer layers, where temperatures are 
around ~ 10"^ — 10'' K. The high values of the creation rates 
in these temperature regimes make H2 increase of a few or- 
ders of magnitude and HD, as well. 

The second test (see Sect. 3.2) is a mean-density cos- 
mological evolution, where the CMB was assumed to be the 
only source of ionizing radiation. We found that (Fig. 6), 
because of the low cosmic background emission, the CMB 

does not affect significantly the abundance evolution, even 
when considering photon propagation at very-high redshift. 

As applications of the numerical methods previously de- 
scribed, we have studied (Sect. 4) the effects of RT on chem- 
ical evolution of primordial gas (e.g. Fig. 7). We performed 
cosmological simulations of structure formation with and 
without ionizing radiation from stellar sources, and checked 
the consequences for molecule formation and destruction. 

We found that the presence of ionizing radiation from 

stars depletes molecular hydrogen up to several orders of 
magnitude (see Fig. 8), by inhibiting the main formation 
paths (the H~-channel and the H^-channel) and by disso- 
ciating it via Lyman- Werner radiation. Our results are con- 
sistent with other studies, e.g. Wise & Abel (2007); Johnson 
et al. (2007); Ahn et al. (2009); Trenti & Stiavelh (2009); 
Whalen et al. (2010); Latif et al. (2011), who studied the 
impacts of the UV background in destroying down the H2 
molecule. In addition, we found that also other molecules 



formed in pristine gas, like e.g. HD, are strongly suppressed 
by radiative feedback (Fig. 9). 

There are large discrepancies in the quantitative as- 
sessments of the impacts of radiative feedback on baryonic 
structure formation in the current literature. In fact, dif- 
ferent works (e.g. Ciardi et al. 2000; Dijkstra et al. 2004) 
show uncertainties of several orders of magnitudes on the 
basic estimates of the masses of the haloes affected by ra^ 
diation. Such determinations could even be highly biased 
by their different "post-processing" approaches. Indeed, in 
order to have large statistics or to circumvent numerical 
complications, these studies do not consider hydrodynam- 
ics self-consistently coupled with RT and non-equilibrium 
chemistry. That could be problematic for getting to reliable 
conclusions, since, as shown by our numerical simulations 
(Sect. 4), RT effects on gas evolution could be quite impor- 
tant and affect star forming regions in a non-negligible way. 
The destruction of early molecules by stellar sources has a 
deep impact on the following star formation processes, be- 
cause it hinders the successive birth of metal-free stars. This 
also implies the need of different viable low-temperature 
coolants, like metals (e.g. Maio et al. 2007), or the presence 
of reionized gas (Yoshida et al. 2007) to sustain star forma- 
tion at later times. The obvious expectation is, then, that the 
popIII star formation rate will drop (even more heavily than 
expected by Maio et al. 2010, 2011; Hasegawa et al. 2009; 
Tornatore et al. 2007; O'Shea & Norman 2008b; Johnson 
& Khochfar 2011). However, further and more detailed in- 
vestigations of high-resolution, numerical simulations are re- 
quired to draw final and definitive conclusions on this topic. 
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